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1. INTRODUCTION 

In this era of emerging global technologies, the global position system (GPS) has advanced 
reasonably in the recent few years. The system allows the user to obtain beneficial and precise information 
regarding the positioning and timing notification respective to a user, with the help of appropriate equipment. 
Since the technology is based on satellite positioning and navigation, it can be accessed globally and 
universally. The messages and ranging information transmitted by the satellite is obtained by the reception of 
the GPS signals, which would be helpful in determining the time and position of the location accordingly. 
From any point on the planet, at any instant of time, at least four satellites can be viewed due to a 24 slot 
arrangement. In order for this to be possible, in the GPS constellation, the satellites are placed in 6 equally 
partitioned planes. Every plane would consist of 4 segments, occupied by the satellites. As the satellites are 
not present in the geosynchronous orbits, the layman would always observe that the geometry of the satellite 
in the visible sky is constantly changing [1], [2]. Many techniques and methods have been studied and 
investigated for identifying outliers in positioning determination using GPS [3]-[6]. Based on least square (LS) 
techniques, many standard receivers have been designed [7]-[9]. All solutions obtained by using the least 
square method would attempt to reduce the sum of squares of the residual errors between the values obtained 
from the estimated models and the values obtained from the measurement. The observations read from the 
equations pertaining to the conventional LS, is presumed to be free of errors, which include both gross and 
systematic errors [10]-[13]. When the measurements of the global positioning system are affected by the 
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outliers, the techniques of LS estimation will result in obtaining an estimation that is bad. So the best method 
to handle outliers is to use a stronger and more robust technique of approach to handle outliers [14]. 
There are a variety of techniques to choose from but not limited to the following such as generalized 
M-estimators, the Danish method, R-estimators, MM-estimators, least absolute values method, M-estimators, 
least median squares, least trimmed squares estimators [3], [4]. 

Moving-horizon estimation (MHE) technique is performing state estimation by using limited and 
most recent information. This technique can be widely used for estimating the state and process monitoring. 
A lot of focus and consideration has been contemplated on MHE due to its rate of success that has been 
achieved for both linear and non-linear systems in predicting dynamic states and parameters [15]-[18]. 
The technique of MHE boasts about using a novel method of collecting only a limited extent of data rather 
than retrieving the entire collection of information from the inception. This technique utilizes only the most 
recent data to find a solution to the optimization technique, where the state estimation is obtained by 
minimizing the least-squares function [15], [16], [19], [20]. In the above-mentioned structure, MHE relies on 
the approach of transforming the problem of estimation using the window of moving fixed size estimation 
into a program that is quadratic. The quadratic program’s size would be determined by the window of fixed 
size estimation. 

GPS signals may contain position errors due to environmental factors or malfunction of subsystems. 
Environmental factors like reflection, attenuation, scattering of signals may contribute to position errors. 
When the sensor or system malfunctions, it may contribute to an error in positioning. This paper addresses 
the position estimation issue of the signal measurements collected from the GPS, which possess outliers. 
Most of the algorithms that we pursue in solving this problem are used in concurrence to [21], and this method 
has in turn been derived from the algorithms that use a methodology of comparison [3], [4]. When the pseudo 
range of a global positioning system signal measurement becomes an outlier, the MHE is used as an estimator to 
determine the measurements. This paper proposes the model based on the algorithms of [1], [21], [22]. 
Then the positions are determined or estimated using the proposed technique. The results obtained from this are 
compared with the results obtained from kalman filtering. The results of the comparison are added to a single 
value of a synthetic outlier, then the position estimation is performed, and the measurements of pseudo range 
are obtained. Using the same concept, the outlier identification algorithm and the MHE are used in parallel 
as in [23], [24], and some other additional measures can be clubbed with the same. The output results 
obtained from the simulation are exclusively designed for the positioning estimation of the measurements 
subjected to outliers and free noise of GPS positioning having the original code [21]. In short, there is no 
new strategy stated here to detect the outliers as such, but the modified moving horizon filter combined with 
the existing GPS parameters would be conjoined together, the result obtained in the process would be 
subjected to comparison with the output obtained from the extended Kalman filter (EKF). 

The paper is organized as follow: in section 2, we introduce the position estimation with a brief 
summary of the error model and sources. Section 3, concerns the extended Kalman filtering model which is 
used for comparison purposes and the proposed moving horizon filter is discussed in section 4. In section 5, 
simulation example is used with results discussed in section 6. Finally, the conclusion is drawn in section 7. 


2. GLOBAL POSITIONING SYSTEM 

Position determination is an issue that has to be solved by the observer, which is implemented with the 
features of receiving and decoding the signals obtained from the GPS. There are 3 dimensions of positioning 
that have to be solved in free space. When the device starts operating autonomously, it is not always presumed 
to be in line with perfect synchronizing with the system time of the satellite. In relevance to the concept 
explained in [1], the non-linear simultaneous equations derived from the system can be used to determine the 
4 variables which would represent measurements obtained from 4 different satellites. 


(X,-—x)* + %-—y)*+ (4, -2)? +b (1) 
pr = Qz- x)? + M-Y} + (Z-z) +b (2) 
ps = VX; =x) + M —y)? + Z3-2)? +b (3) 
(X,— x)? + Y4- y)? + (Z4 —2)* +b (4) 
Where p;: noiseless pseudo-range, (X; Y; Z;)': coordinates of satellite i (J = 1,...,4), (x y zZ)": user 


coordinates, and b: receiver clock. To get the four known parameters, on the (1)-(4) are linearized about an 
approximate observer location Xo, = (Xo Yo Zo)! fori = 1,..., 4. 
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RE E o E E 
ax (Xj-X9)*+ (Vi-Yo)?+ (Zi-Zo)” (5) 
OE a A SE 
dy (Xj-X0)*+ (Yi-Vo)2+ (Zi-Zo)? (6) 
ON oa ye E E 
dz I(Xi-X0)*+ Vi Vo)? + Zi-Zo) (7) 


The resulting measurement vector equation for pseudo-range as the observable is given in (8) as a noiseless 
model; where f;(x,,): predicted pseudo-range, (Ax Ay Az)”: difference vector, x,, is the nominal point 
linearization based on (xo Yo Zo)! and predicted receiver time. 


db. Ob. ddr 4 
Ox Oy 02 
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X _| ax y az Ay 
Pi — Pi(Xor) = av, avs avs 4 || az (8) 
Ox Oy Oz b 
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P41 
_ | P2 
P= | ps (9) 
P4 
P1(%or) 
zs P2(Xor) 
(Xe) Se 2 10 
Be ae (10) 
Pa(Xor) 


The noiseless pseudo range can be written as 


pi = V(Xi-— x)? + Yi- y)? + (Ziz)? +b (11) 


2.1. GPS error model and sources 

Errors in measurements are classified into noise and bias. Although they both seem the same, noise 
is a type of error that will change rapidly and returns to a null value within a quick interval. Bias is a type of 
error that may exist over a period of time. The errors in GPS may occur due to various factors. 
The navigation message that is broadcasted by the satellite may contain some errors in the values of the 
parameters can cause errors. When the signal is traveling from the satellite to the receiver, the signal 
propagation medium varies based on some uncertainties, which affects the signal’s traveling time which may 
contribute to error. Another source of error may be caused due to the receiver’s noise [1], [2], [22], [25]-[27]. 
Some of the significant errors likely to affect a GPS model, like control system error, signal propagation 
error, and measurement errors, are discussed below. 


2.1.1. Control system error 

The cause of this type of error is based on the measurements of the satellite clock and the ephemeris. 
They are calculated at the GPS monitoring station in the control segment. With the utilization of a Kalman 
filter, the satellite position, satellite velocity, clock bias, frequency bias, and many other such parameters are 
obtained. Along the three orthogonal directions, the decomposition of the ephemeris error is carried out. 
Root sum square value of the clock error and line of sight of ephemeris error contributes to the range error of 
errors in clock and ephemeris [22], [28], [29]. 


2.1.2. Signal propagation modeling error 

These categories of errors are mainly caused due to the changes in carrier phase and code 
measurement values [1], [22]. Which can briefly: 1) signal refraction, dispersive medium, and wave 
propagation: when the signal of the GPS travels from the satellite to the receiver, the medium of propagation 
affects the signal properties. Since the signals pass through various layers of the atmosphere, it passes 
through the ionosphere which consists of charged particles, and it passes through the troposphere which is the 
neutral portion of the atmosphere. This non-uniformity in the medium of transmission would change the 


TELKOMNIKA Telecommun Comput El Control, Vol. 20, No. 2, April 2022: 426-436 


TELKOMNIKA Telecommun Comput El Control o 429 


refractive indices, which may alter the speed of the signal, thus causing errors; and 11) ionospheric and 
tropospheric delay: The ionosphere tends to be generally stable, but it may tend to exhibit certain fluctuations 
at areas that are near the poles and the equators. The electron densities tend to vary and change as a result of 
certain solar storms or flares that can cause errors in signal tracking, as carrier phase and amplitude change. 
In the troposphere, a delay might be experienced by the signal as a result of the refractive index of air mass. 
The density of the water vapor changes in accordance with the weather condition, in such conditions the 
modeling becomes difficult. Also, a dry atmosphere causes a delay [30]-[32]. 


2.1.3. Measurement error 

These types of errors are caused by the power of the signal, the structure of code, and the design of 
the antenna and receiver. Multipath error, receiver noise error, errors due to measurement models, and user 
ranges can affect the signal and cause outliers during the process of estimation. Such types of errors are 
briefed as follow: 1) receiver noise and multipath: when a random noise smudges the signal that is responsible 
for the precision of the code and carrier phase measurement, it is referred to as receiver noise. Many factors 
may lead to this type of error; some of the factors are antennas or amplifier noise, cable/receiver noise, noise 
due to signal quantization. If there is no interfering noise present, then the receiver waveform is a sum of the 
received GPS signal and the fluctuating noise that is random [1], [22], [33]. Multipath causes a change in 
measurement of phase due to an interfering signal. When the signal approaches the antenna through more than 
2 or more paths or ways, the antenna will receive multiple signals, which includes the direct line of sight signal 
and the various reflected signals. Usually, the multipath error is quasi sinusoidal [33]-[35]; and 11) measurement 
error model and user range error: the source of the error model is the receiver noise and multipath, which 
depends on the elevation angle. When the satellite gets to the lower levels of the atmosphere, the power signal 
reduces, and hence the multipath error tends to increase, which can be controlled with high-end receivers. 
User range error is the square root of the sum of squares of the range error of control segment, atmospheric 
propagation, and receiver noise and multipath [30], [36]. 


3. ESTIMATION STRATEGY AND OUTLIER DETECTION 

Estimating the state variables when a system is contaminated by outliers can be treated using 
various approaches, including MHE and Kalman filtering. As it is well known, Kalman filter is the best 
estimator in terms of minimizing the quadratic estimation error with white Gaussian disturbances. Extended 


Kalman filter (EKF), which is the non-linear version of the Kalman filter and the model can be described 
in (12.a) and (12.b). 


Xt = f (Xt-1) + Wee (12.a) 
Z, = hX) +», (12.b) 


Where X, is the state, Z, is the measurement, f(X;_,) and h(X;) are the state and measurement functions, 
respectively, w;_,1S system noise with covariance Q, and v, is measurement noise with covariance R. 
EKF prediction and update are performed in two stages; 
— Prediction 

State estimate: Šit- = f Xe-ae-1) 

Error covariance: Pyt- = Fey Pr—aje-1 Fea + Qe- 
— Update 

Kalman gain: K; = (Pyye-1 H7 (Ae Peje-1 H7) + RD 

State estimate: Rijt = Šit- + K (z — h(Rrt-1)) 

Covariance estimate: Pae = U — KH) Pejt-1 

Extended Kalman filter approximate non-linear function, whereby linearized kalman filter can be used. 
Figure 1 shows the Kalman filter setup. MHE is a topic that gathered a high degree of contemplation riding on the 
success model which has been achieved with predictive control. Based on the information proposed by Jazwinski 
in his pioneering work [37], the basic idea used in MHE for dynamic systems is to perform a state estimation 
utilizing the limited quantity of most current or recent information. The most important benefit of using the MHE, 
is its capacity to be accountable for constraints that are present on the state variable [15], [16], [19], [20], [38]-[40]. 
Outliers prevent the ensuring of a guaranteed performance of an estimator by adding a particular kind of 
uncertainty in general [41]. In the process of designing filters for any uncertain system, the key requirement 
to be remembered is the robustness of the system. 
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COPY OF 
DISCRETE SYSTEM 


Figure 1. Kalman filter set-up 


4. MODIFIED MOVING-HORIZON FILTER FOR POSITIONING ESTIMATION WITH 
PSEUDO-MEASURE SUBJECT TO OUTLIER 
Alessandri et al. [23], [24], have developed a novel method for states variable estimation with 
outlier affecting the measurements equations by minimizing least-squares cost function, where the strategy is 
to leave out in turn all the measurement that can be affected by the outlier. The design will be used in our 
simulation considering the non-linear approach, which can be presented for the case in which measurement is 
corrupted by outlier as follow. 


a Zz = 1 Z 
Cra a H4 \Z:-v.e— Xen ll? i=) : EN: OF MERS (13) 
N 


i+t-N+k-1 
For k = 1,2,...,N +1 


where u > 0, the cost in (13) is to be minimized together with the constraints £i41,t = f (it ui) 
for it=t-—N,..t—1. Practically, at each tme t=N,N+1,... we need to solve 
MİNG s.t, 2:44 =f (finu) JE(®) for k = 0,1,...,N +1 and % € R”. The estimates of x,_y,...X; at time t is 
represented as £t-y to ft-N+1,t =- tt - By comparing the optimal cost of both cases (no outlier case and 
outlier affecting the measurement case), the best cost is associated with estimate, and the approach would be 
presented as: 


: : Kra 
PUN rao tnd PUN Sst tiit S rad) Je) (14) 
at each time t = N,N +1, ... 


Such algorithm considering the output £+-y, t, can be described as shown in Table 1. The strategy is 
implemented for simulation purposes following the algorithm in Table 1 and in line with [23], [24], 
the flowchart of the basic proposed moving-horizon filter (MHF) set up at time t is shown in Figure 2. 


Table 1. Outlier MHF algorithm 


No. Description  _—_—_— — 
1 fork =0toN +1 compute 2E y: 


Ririt = f (iz, ui) holds for argmin Jf (Z;) 


2 end for 
3 Choose kj € argmin JE (2,-n,t) 
4 * 


. x 2. SIG 
Update estimate X,_yj2 = Xp yit 


As shown in Figure 2, the algorithm works on minimizing a set of cost functions (LS), by leaving out all 
possible outliers within the bath, where at each time instant the lowest cost is determined and hence pushed 
forward and the process is repeated for all batches (t = N, N+1,...). 


5. METHODOLOGY 
In this section we address the simulation of position estimation when pseudo-measure 1s 
contaminated by outlier. The MHF and EKF proposed in section 3 and 4, respectively, are simulated for the 
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case where no pseudo-measures are contaminated by outlier, and for the case of outlier occurrence (1.e. single 
outlier at fixed position). The principle followed here is in line with [1], [2], [21], and considering the 


previous equations. The system is modeled in (15). 


pi =~ptbt+y; (15) 
Where ¢ọ = 4/ (X; —x)* + (j-—y)*+ (Z,-—2)?, clock bias b, and {v;} is unknown pseudo-range 
measurement noise which is zero-mean Gaussian except in case of outliers occurrence, which unpredictable 
and has an unknown and very large covariance. The state, can be defined as in (16). 


X=(X ù% Yy Z Ub qd) (16) 


Where position is represented in x, y, Z coordinate, Vy, Vy, and v, are the velocity, and d is the clock drift. 
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Figure 2. MHF implementation at time t 


5.1. Case study 
A case of the estimation of position combined with a pseudorange measurement is taken as an outlier, 


is considered here. The original sampling data are taken from [42], and the noiseless model is from [21]. Let the 
covariance matrix be defined as: 


Q = (Qx Qy Qz Qp) (17) 
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p3 T7 
Sox + Sox 3 Si Eg 
Qy = X (18) 
Sox > SoxT 
T? T42 
SeT +S, — S,— 
f g g 
Qs = a : (19) 
Sg oe Sof 


Where Spx is the power spectral density of speed noise, Sp is the power spectral density of clock bias noise, 
and Sp is the power spectral density of frequency drift noise. The initial values of X and f(.) are taken 
according to [1], [21], and the setting of the measurement noise covariance R is refereed to [43]. The setting 
of MHF is the same as presented in [23] with additive running criteria. We simulate the pseudo-range and 
satellite position of a GPS receiver at fixed location for a period of 25 seconds. Synthetic outlier has been 
added to a pseudo-range measurement at receiving time t = 9 s. For the comparison with EKF, the estimates 
of the moving-horizon filter obtained by solving £+—y + with the initial setup briefed in Table 2. 


Table 2. MHF estimates set-up 
Ci eseription > 


mn — akt) 
Solve X:_ne = Xp 


ki € argmin Jf (8-2) fork = 0,1,..,N+1 


Estimates obtained with 
N = 4 window size 
u = 0.6 Tuning parameter 


For 25 simulation run (i.e T = 25 steps) 


6. RESULTS AND DISCUSSION 

The performance of the proposed modified MHF with tuning parameters described in Table 2 is 
compared with EKF and LS techniques. The positioning estimation where no measurements are affected by 
outliers (i.e simulating the original sample with zero mean Gaussian noise) is shown in Figure 3(a) 
for a period of 25 seconds. The positioning error is shown in Figure 3(b) shows that the maximum error 
recorded for MHF occurred at first 3 readings (due to initial value) with median values 1.45, -0.42, and -2.7 
for coordinate (x, y, Z) respectively for 25 readings. The case in which the pseudorange measurements are 
contaminated by outliers is shown in Figure 4 for a case of a synthetic outlier with a random value and fixed 
position at t= 9 s. The original value is kept in the plotting for comparison reasons, where the arrow indicates 
the synthetic outlier. The simulation results show that the proposed MHF is more robust to outliers than EKF 
where the MHF has not been so affected by outlier, still the effect of outliers is observed for the 
measurements following the time step of the outlier occurrence, which can be improved with more simulation 
runs and careful choice of the tuning parameter. As shown in Figure 4(a) position estimation [x, y, z] and 
Figure 4(b) the maximum relative positioning errors were captured for LS and EKF at the step of outlier 
occurrence and t + 3 following steps. The maximum relative errors for each coordinate are shown in Table 3. 


Table 3. Maximum relative error 
Step Coordinate Filter Relative error 


11 X LS 93.31 
O loutlier Y EKF 109.37 
=l Z -ise e 


Compared with the case of no outliers, the EKF is quite sensitive to outlier value, the highest errors 
captured at time instant of outlier occurrence, and the following time steps where outlier effect is observed. 
From the computational time point of view, the total time elapsed is around 60 seconds for 25 simulation 
runs (for all estimators) with mean computational time (MCT) of EKF around 0.008 s, versus 0.028 s for 
MHF (the computational time of MHF turns to be demanding comparing to EKF, but the performance turns 
to be similar over different simulations opposite to EKF which is affected by outlier position, outlier value, 
and number of runs), such results are varied with window size (N). 
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Figure 3. The case of no outlier: (a) position estimation [x, y, z] and (b) relative positioning error [x, y, z] 
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Figure 4. The case of outlier [t = 9]: (a) position estimation [x, y, z] and (b) relative positioning error [x, y, z] 


7. CONCLUSION 

This paper proposes applying the concept of moving horizon estimation to the problem of state 
estimation in which pseudo range measures are affected by outliers. With performed simulation, the results have 
proved that EKF is turning to be sensitive to outliers opposite to the proposed MHF. The ability of the proposed 
filter to leave out a measure that is possibly contaminated by outlier is proved by the simulation mean and 
provided algorithm. The proposed estimator can be easily tuned with the correct selection of tuning parameter, 
however, the computational time of such estimator is still challenging, which will be studied in future work as 
well as the effort in reducing the computational time. 
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